#command line arguments
#args <- commandArgs(TRUE)

#reading the beagle utility output file of the window based
#data<-read.table("mb1_beagle_util",sep="\t")  #change paramteres
#data<-read.table(as.character(args[1]),sep="\t")

data<-read.table("halfmb_beagle_util",sep="\t")
colnames(data)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
distance<-read.table("../1000genome/genome1000_4mb1mbedge/1/markers.txt",sep=" ")
colnames(distance)<-c("marker_identifier","dist","allele1","allele2")
merge1<-merge(data,distance,by="marker_identifier")

data<-c()
data<-read.table("halmb_s1_beagle_util",sep="\t")
colnames(data)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
merge6<-merge(data,distance,by="marker_identifier")

data<-c()
data<-read.table("halfmb_s2_beagle_util",sep="\t")
colnames(data)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
merge7<-merge(data,distance,by="marker_identifier")

#reading the beagle utility output file of the whole chr imputation
data<-c()
data<-read.table("mb4_beagle_util",sep="\t")
colnames(data)<-c("marker_identifier","minor_allele","minor_allele_frequency","allelic_r_squared","dosage_r_squared","HWE_dosage_r_squared","accuracy","missing_score")
merge3<-merge(data,distance,by="marker_identifier")


#file names
window<-"halMB_SEED1_SEED2"  #change paramteres
#window<-as.character(args[2])
allelicr2_name<-paste("whole_default_vs_",window,"_allelic_r_squared.pdf",sep="")
dosager2_name<-paste("whole_default_vs_",window,"_dosage_r_squared.pdf",sep="")
HWE_dosage_r_squared_name<-paste("whole_default_vs_",window,"_HWE_dosage_r_squared.pdf",sep="")

#parameters for this run
window_start<-c(16050035,18050036,20050036,22050036,24050036,26050036,28050036,30050036,32050036,34050036,36050036,38050036,40050036,42050036,44050036,46050036,48050036,50050036)
window_stop<-c(18050036,20050036,22050036,24050036,26050036,28050036,30050036,32050036,34050036,36050036,38050036,40050036,42050036,44050036,46050036,48050036,50050036,51244207)
actual_window_start<-c(16050035,17550036,19550036,21550036,23550036,25550036,27550036,29550036,31550036,33550036,35550036,37550036,39550036,41550036,43550036,45550036,47550036,49550036)  #change paramteres
actual_window_stop<-c(18550036,20550036,22550036,24550036,26550036,28550036,30550036,32550036,34550036,36550036,38550036,40550036,42550036,44550036,46550036,48550036,50550036,51244207)  #change paramteres
#para_start <- as.character(args[4])
#para_stop <- as.character(args[5])
#actual_window_start<-as.numeric(unlist(strsplit(para_start,',')))
#actual_window_stop<-as.numeric(unlist(strsplit(para_stop,',')))





num_windows<-length(window_start)
edge<-500000      #change paramteres
merge2<-merge1[which(merge1$dist >=window_start[1]+edge & merge1$dist <= window_stop[num_windows]-edge),]
merge1<-c()
merge4<-merge3[which(merge3$dist >=window_start[1]+edge & merge3$dist <= window_stop[num_windows]-edge),]
merge3<-c()
merge8<-merge6[which(merge6$dist >=window_start[1]+edge & merge6$dist <= window_stop[num_windows]-edge),]
merge6<-c()
merge9<-merge7[which(merge7$dist >=window_start[1]+edge & merge7$dist <= window_stop[num_windows]-edge),]
merge7<-c()



#creating the indicator for different windows of window based
merge2$ind<-0
i<-1
while(i<num_windows+1)
{
	merge2$ind[which(merge2$dist >=window_start[i] & merge2$dist <=window_stop[i] )]<-i	
	i<-i+1
}

#creating the indicator for different windows of whole chr
merge4$ind<-0
i<-1
while(i<num_windows+1)
{
	merge4$ind[which(merge4$dist >=window_start[i] & merge4$dist <=window_stop[i] )]<-i	
	i<-i+1
}

#creating the indicator for different windows of whole chr
merge8$ind<-0
i<-1
while(i<num_windows+1)
{
	merge8$ind[which(merge8$dist >=window_start[i] & merge8$dist <=window_stop[i] )]<-i	
	i<-i+1
}


#creating the indicator for different windows of whole chr
merge9$ind<-0
i<-1
while(i<num_windows+1)
{
	merge9$ind[which(merge9$dist >=window_start[i] & merge9$dist <=window_stop[i] )]<-i	
	i<-i+1
}




#calculating distance from edges for the window based
nc<-nrow(merge2)+1
i<-1
distance<-c()
while(i<nc)
{
#	print(i)
	x<-abs(merge2$dist[i]-actual_window_start[merge2$ind[i]])
	y<-abs(merge2$dist[i]-actual_window_stop[merge2$ind[i]])
	if(x <y)
	{
		distance[i]<-x
	}else{
 		distance[i]<-y
	}
	i<-i+1
}
merge2$distance<-distance
distance<-c()

#calculating distance from edges for the whole chr based seed1
nc<-nrow(merge8)+1
i<-1
distance<-c()
while(i<nc)
{
#	print(i)
	x<-abs(merge8$dist[i]-actual_window_start[merge8$ind[i]])
	y<-abs(merge8$dist[i]-actual_window_stop[merge8$ind[i]])
	if(x <y)
	{
		distance[i]<-x
	}else{
 		distance[i]<-y
	}
	i<-i+1
}
merge8$distance<-distance
distance<-c()

#calculating distance from edges for the whole chr based seed2
nc<-nrow(merge9)+1
i<-1
distance<-c()
while(i<nc)
{
#	print(i)
	x<-abs(merge9$dist[i]-actual_window_start[merge9$ind[i]])
	y<-abs(merge9$dist[i]-actual_window_stop[merge9$ind[i]])
	if(x <y)
	{
		distance[i]<-x
	}else{
 		distance[i]<-y
	}
	i<-i+1
}
merge9$distance<-distance
distance<-c()

#calculating distance from edges for the whole chr based
nc<-nrow(merge4)+1
i<-1
distance<-c()
while(i<nc)
{
#	print(i)
	x<-abs(merge4$dist[i]-actual_window_start[merge4$ind[i]])
	y<-abs(merge4$dist[i]-actual_window_stop[merge4$ind[i]])
	if(x <y)
	{
		distance[i]<-x
	}else{
 		distance[i]<-y
	}
	i<-i+1
}
merge4$distance<-distance
distance<-c()

#plotting allelic r2
pdf(allelicr2_name,width =50,height =50)
plot(merge2$distance, merge2$allelic_r_squared, xlab="Distance from edge", ylab="allelic_r_square", pch='.')
points(merge4$distance, merge4$allelic_r_squared, pch='.')
points(merge8$distance, merge8$allelic_r_squared, pch='.')
points(merge9$distance, merge9$allelic_r_squared, pch='.')
lines(lowess(merge2$distance,merge2$allelic_r_squared), col="purple") 
lines(lowess(merge4$distance,merge4$allelic_r_squared), col="green")
lines(lowess(merge2$distance,merge8$allelic_r_squared), col="red") 
lines(lowess(merge4$distance,merge9$allelic_r_squared), col="blue")
legend("top", horiz=T,col=c("purple","green","red","blue"), lty=c(1,1,1,1), c("Half MB edge 2MB window lowess curve","4mb window 2mb edge lowess curve","Half MB edge 2MB window seed 1 lowess curve","Half MB edge 2MB window seed 2 lowess curve"), pch=c('.','.',-1,-1),cex=3)
dev.off()

#plotting dosage r2
pdf(dosager2_name,width =50,height =50)
plot(merge2$distance, merge2$dosage_r_squared, xlab="Distance from edge", ylab="allelic_r_square", pch='.')
points(merge4$distance, merge4$dosage_r_squared, pch='.')
points(merge8$distance, merge8$dosage_r_squared, pch='.')
points(merge9$distance, merge9$dosage_r_squared, pch='.')
lines(lowess(merge2$distance,merge2$dosage_r_squared), col="purple") 
lines(lowess(merge4$distance,merge4$dosage_r_squared), col="green")
lines(lowess(merge2$distance,merge8$dosage_r_squared), col="red") 
lines(lowess(merge4$distance,merge9$dosage_r_squared), col="blue")
legend("top", horiz=T,col=c("purple","green","red","blue"), lty=c(1,1,1,1), c("Half MB edge 2MB window lowess curve","4mb window 2mb edge lowess curve","Half MB edge 2MB window seed 1 lowess curve","Half MB edge 2MB window seed 2 lowess curve"), pch=c('.','.',-1,-1),cex=3)
dev.off()

#plotting HWE dosage r2
pdf(HWE_dosage_r_squared_name,width =50,height =50)
plot(merge2$distance, merge2$HWE_dosage_r_squared, xlab="Distance from edge", ylab="allelic_r_square", pch='.')
points(merge4$distance, merge4$HWE_dosage_r_squared, pch='.')
points(merge8$distance, merge8$HWE_dosage_r_squared, pch='.')
points(merge9$distance, merge9$HWE_dosage_r_squared, pch='.')
lines(lowess(merge2$distance,merge2$HWE_dosage_r_squared), col="purple") 
lines(lowess(merge4$distance,merge4$HWE_dosage_r_squared), col="green")
lines(lowess(merge2$distance,merge8$HWE_dosage_r_squared), col="red") 
lines(lowess(merge4$distance,merge9$HWE_dosage_r_squared), col="blue")
legend("top", horiz=T,col=c("purple","green","red","blue"), lty=c(1,1,1,1), c("Half MB edge 2MB window lowess curve","4mb window 2mb edge lowess curve","Half MB edge 2MB window seed 1 lowess curve","Half MB edge 2MB window seed 2 lowess curve"), pch=c('.','.',-1,-1),cex=3)
dev.off()


